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Unlike the bewildering situation in the 77 * —>• n form factor, a widespread view is that pertur¬ 
bative QCD can decently account for the recent BaBar measurement of 77 * —>• r/ c transition form 
factor. The next-to-next-to-leading order (NNLO) perturbative correction to the 77 * -A r/ c ^ form 
factor, is investigated in the NRQCD factorization framework for the first time. As a byproduct, 
we obtain by far the most precise order-a 2 NRQCD matching coefficient for the r/ Ct i -A 77 process. 
After including the substantial negative order-a 2 correction, the good agreement between NRQCD 
prediction and the measured 77 * -A rj c form factor is completely ruined over a wide range of momen¬ 
tum transfer squared. This eminent discrepancy casts some doubts on the applicability of NRQCD 
approach to hard exclusive reactions involving charmonium. 

PACS numbers: 12.38.Bx, 13.60.Le, 14-40-Pq 


The advent of high-luminosity e + e _ colliders affords an 
invaluable opportunity to study hard exclusive processes 
involving light mesons such as TT,r),r]'. New measure¬ 
ments constantly emerge, some of which greatly challenge 
our understanding toward perturbative QCD (pQCD). 
In 2009, BaBar collaboration reported their measure¬ 
ment of the 77* -A 7T° form factor for 4 GeV 2 < 
Q 2 < 40 GeV 2 Q. Surprisingly, the measured value 
of <3 2 |-F(Q 2 )| increases with Q 2 , exceeding the asymp¬ 
tote v^/tt = 0.185 GeV predicted in pQCD s after 
Q 2 > 10 GeV 2 . This triggered renewed theoretical inter¬ 
est in scrutinizing the leading-twist pion light-cone dis¬ 
tribution amplitude (LCDA), and it looks ironical that 
to date our knowledge of pion is still rather inadequate. 

Apart from hard exclusive reactions involving light 
hadrons, exclusive charmonium production also emerge 
as a new frontier of pQCD in recent years [§]. This topic 
was initiated by the double-charmonium production pro¬ 
cess e + e _ -A J/ijj+r] c first observed by Belle in 2002 Q. 
Disquieting discrepancy between the data and the lead¬ 
ing order (LO) pQCD predictions HH3 has stimulated a 
great amount of theoretical activity. 

In this work, we are interested in examining the 
simplest exclusive charmonium production process: the 
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7*7 Vc transition form factor, which is defined through 
(^c(p)|<?emI7(M)) = ie 2 e fJ ’ vprT £ u q p k (7 F(Q 2 ), (1) 


where is the electromagnetic current, and the mo¬ 
mentum transfer squared q 2 = (p — k ) 2 = — Q 2 < 0. 

In 2010 BaBar collaboration presented very precise 
measurements for the 7*7 -A r]c transition form factor in 
the range 2 GeV 2 < Q 2 < 50 GeV 2 @|. The measured 
F{Q 2 ) can be well described by a simple monopole fit: 

= TT^Ta’ (2) 


with the pole parameter A = 8.5 ± 0.6 ± 0.7 GeV 2 . In¬ 
terestingly, a variety of phenomenological studies, i.e., 
k± factorization 0, H3, lattice QCD [u|, J/ip- pole- 
dominance 0, QCD sum rules (l2 |. light-front quark 
model [lil , all yield predictions compatible with the data, 
at least in the small Q 2 range. 

However, we caution that the satisfactory agreement 
between the existing predictions and data should not be 
interpreted as a triumph of pQCD. On the contrary, as 
we will see later, understanding the 7*777 c form factor 
actually poses an outstanding challenge to pQCD. 

Both light-cone approach 0 and nonrelativisic QCD 
(NRQCD) factorization [l3| are well-founded pQCD 
methods to tackle hard exclusive reactions involving 
charmonium. The bulk of r] c events produced in two- 
photon fusion concentrates in the range Q 2 ~ m 2 , where 
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light-cone approach is obviously inappropriate. Never¬ 
theless, since the charm quark propagator carries typi¬ 
cal virtuality of order m 2 3> Aq CD even in the Q ► 0 
limit, the NRQCD approach, which is based on effective 
field theory and fully exploits the nonrelativistic nature 
of quarkonium, remains as a valid tool. As Q m, both 
light-cone and NRQCD approaches can apply, but lack 
of precise knowledge of charmonium LCDA severely hin¬ 
ders the predictive power of the former approach. In con¬ 
trast, NRQCD approach is much more predictive, since 
the only nonperturbative inputs are a few NRQCD ma¬ 
trix elements that are strongly constrained by the mea¬ 
sured charmonium annihilation decay rates. The consen¬ 
sus now is that the NRQCD approachgenerally yields 
less ad-hoc, more reliable predictions |3|, ll5(, in partic¬ 
ular amenable to systematically incorporating perturba¬ 
tive and relativistic corrections. 

The aim of this work is to investigate higher-order per¬ 
turbative corrections to the 7*7 —> rj c form factor in the 
NRQCD framework. This is motivated by the fact that 
QCD radiative corrections in exclusive charmonium pro¬ 
cesses are often substantial, exemplified by the sizable 
next-to-leading order (NLO) QCD correction to e + e _ —»■ 
J/tp + ri c [16H18I |. and substantial next-to-next-to-leading 
order (NNLO) perturbative corrections to a few 5-wave 
quarkonium decay channels, e.q., J/ip — > e + e~ jH, 0], 
Vc — t 77 0, and B c —> Iv 22, [23|. To the best of our 
knowledge, this work constitutes the first comprehensive 
NNLO analysis for exclusive quarkonium production pro¬ 
cess. 

According to the NRQCD factorization 0 , the form 
factor at LO in v can be expressed as 


F(Q 2 ) = C(Q,m,n Rl fj, A ) 


(Vc\^x(VA)\0) 


m 


+ 0(v 2 ). (3) 


C(Q,m, hr, n A ) represents the NRQCD short-distance 
coefficient, where m, hr, /jl A denote the charm quark pole 
mass, the renormalization scale, and the NRQCD factor¬ 
ization scale, respectively. In phenomenological studies, 
this nonpertubative matrix element is often substituted 
by y/2N c 'ip Vc (0) (N c = 3 is the number of color), where 
ip rit , (0) is the Schrodinger wave function at the origin for 
rj c in quark potential model. In NRQCD, this matrix 
element becomes a scale-dependent quantity. 

Since the vacuum-to-^ c matrix element cancels in the 
ratio F(Q 2 )/F( 0), the NRQCD prediction to ([2]) is free 
from any adjustable nonperturbative inputs. In some 
sense, this normalized r] c form factor bears some similar¬ 
ity with the ratio T[J/tp e + er]/T[ r q c —> 77 ], which has 
been carefully studied in [21, 241. However, for the lat¬ 
ter ratio, the NRQCD matrix elements of J/ip and r/ c do 
not exactly cancel due to the presence of the spin-spin 
interaction, the estimate of which requires the nonper¬ 
turbative ansatz of modelling the charmonium. In this 
respect, confronting the BaBar data of the r/ c form fac¬ 
tor over a wide range of Q 2 would pose a unique and more 


nontrivial test against the prediction from NRQCD fac¬ 
torization. 

Owing to the fact Q,m Aqcd, the short-distance 
coefficient can be reliably computed in perturbation the¬ 
ory. Through 0(a 2 ), it can be organized as 


C(Q,m,n R ,n A ) =C ( 0 ) (Q,m)jl + C' F 
a? fAi.in /4 Cf/(i) (t) - n 2 C F 


4 Q 2 + m 2 


x In — + /^(t) 
m 


+ 0(aM, 


as ^) /(i)( r ) 
7 r 



( 4 ) 


where r = Q 2 /m 2 . C F = f and Ca = 3 denote the 
51/(3) color factors. 

(7(°) signifies the LO NRQCD matching coefficient: 


C* 0 >(Q,m) 


4eg 

Q 2 + 4to 2 ’ 


( 5 ) 


where e c = | (e& = — |). As expected, C oc 1 /Q 2 
asymptotically. The LO NRQCD prediction to F(Q 2 ) 
coincides with the monopole form in m, with the re¬ 
spective pole parameter A = 4m 2 ~ 9 GeV 2 . 

The 0(a s ) contribution is encapsuled in /A)(r): 


f W (r) 


7t 2 (3 — t) 20 + 9r r (8 + 3r) 4 + r 

6(4 + r) ~~ 4(2 + r) ~ 4(2 + r ) 2 2 



T 

2(4+ r) 




( 6 ) 


which can be translated from Eq. (A3) of j25| . 

In Eq. @, /3o = ^ C A — | T f (til + n F ) is the one- 
loop /3-function coefficient, n F = 1, and til labels the 
number of light quark flavors (til = 3 for rj c , 4 for 775 ), 
and T F = The occurrence of the /3 0 In term in 
P|) is dictated by renormalization-group invariance. The 
occurrence of the factorization scale p, A reflects that the 
NRQCD pseudoscalar density tp'x has a non-vanishing 
anomalous dimension starting at order-a 2 , as first dis¬ 
covered in 0 - The nontrivial task is to decipher the 
function /^(r). 

We start with computing the on-shell quark amplitude 


for 7*7 —> cc^S^) through order a 2 , using the covariant 
trace technique to facilitate the projection of the cc pair 
onto the intended quantum number. Prior to performing 
the loop integration, we neglect the relative momentum 
between c and c, which amounts to directly extracting 
the short-distance coefficients at v° accuracy H, 20l |. 

We briefly outline the calculation. We employ the 
package FeynArts [26] to generate corresponding Feyn¬ 
man diagrams and amplitudes through order-a 2 in Feyn¬ 
man gauge. There are 108 regular 2-loop diagrams and 
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FIG. 1: Sample Feynman diagrams for 7*7 —► cc( 1 S'q 1 ^). 


12 “light-by-light” scattering diagrams, two of which are 
sketched in Fig. [Q The latter gauge-invariant subset is 
UV- and IR-finite. Dimensional regularization (DR) is 
utilized to regularize both UV and IR divergences. We 
then use FeynCalc/FormLink [27], HU to carry out the 
trace over Dirac and color matrices. We follow [29] on 
the prescription of 75 in DR. The packages Apart (3d ] 
and FIRE 31i] are utilized to conduct partial fraction 
and integration-by-parts reduction. Finally, we end up 
with the 174 Master Integrals (MI). For a handful of 


simpler Mis, we employ the a parameters m as well 


as the Mellin-Barnes tools AMBRE 


1 /MB 


341 ] to infer 


the (semi)-analytic expressions, while for the remaining 
Mis, we combine FIESTA CubPack [35, 36j to carry out 
sector decomposition and subsequent numerical integra¬ 
tions with quadruple precision. 

In implementing the quark wave function and mass 
renormalization, we take the order-a 2 expressions of Z 2 


and Z m from [37, |38|. The strong coupling constant is 
renormalized to one-loop order under MS scheme. All 
the UV divergences are removed after the renormaliza¬ 
tion procedure. At this stage, the NNLO amplitude still 
contains a single IR pole with the very coefficient as an¬ 
ticipated from ([4]). We factorize this IR pole into rj c - to- 
vacuum NRQCD matrix element according to MS pre¬ 
scription, with In pa manifested as the remnant. Ulti¬ 
mately, we are able to identify each individual term as 
specified in (|4]), with high numerical precision. This is 
the first calculation that endorses the validity of NRQCD 
factorization at the NNLO level for exclusive quarkonium 
production process. 


It is convenient to divide the NNLO contribution 
/ ( 2 ) (r) in ffj into the regular piece and the “light-by¬ 
light” piece, / ( 2 ) (t) = /reg(r) + //^(r). The former is 
real-valued, and the latter is complex-valued, whose sig¬ 
nificance and computational cost explode with increasing 
r. /reg (t), as plotted in Fig. [2j turns to be deeply nega¬ 
tive, signaling a substantial negative NNLO correction to 
the form factor. Despite being a rather flat function of r, 
the asymptote of //eg (t) appears to be numerically con¬ 
sistent with the 2 -loop leading collinear logarithm In 2 r 
predicted in @. Some benchmark values of //<?g (r) and 
fibi ( T ) are listed in Table [T] In the experimentally ac¬ 
cessible range of Q 2 , the “light-by-light” contribution for 
r] c appears insignificant relative to the regular one. Nev¬ 
ertheless, the “light-by-light” piece turns out to be non- 
negligible even at relatively low r for 


T 

1 

5 

10 

25 

50 

f (2) 

J reg 

-59.420(6) 

-61.242(6) 

-61.721(7) 

-61.843(8) 

-61.553(8) 

-f (2) 

Jlbl 

0.49(1) 

-0.48(1) 

-1.10(1) 

-2.13(1) 

-3.07(1) 

—0.65(l)i 

—0.72(l)i 

—0.71(l)i 

—0.69(l)i 

-0.68(1)* 

r(2) 

J reg 

-59.636(6) 

-61.278(6) 

-61.716(7) 

-61.864(8) 

-61.668(8) 

A 2 ) 

Jlbl 

0.79(1) 

-5.61(1) 

-9.45(1) 

-15.32(1) 

-20.26(1) 

— 12.45(1)* 

—13.55(1)* 

—13.83(1)* 

-14.03(1)* 

-14.10(1)* 


TABLE I: frel (r) and (r) at some typical values of r. The 
first two rows for r/ c and the last two for r/b- 


The t] c form factor at zero momentum transfer squared 
reads 

F <°) = A/, Wx(m)| 0 >{ 1 + (T - i) 


at 


Cf 


A _ in A -2 


2/4 m 2 


- 7r CV C f + 


Ca \ , Ma 


In ■ 


+/r ( ei(0) + ft b7(0) 




(7) 


where 



/reg (0) = -21.10789797(4)C'| - 4.79298000(3)C F CU 

W + 5 ln2 + i c( 3 ) ^S) CFTFni 

+0.223672013(2)U F T F n H , ( 8 ) 


and 


/$(0) = (V73128459 + in - 0 )c F T F £ 

0.64696557 + 2.07357556 i]c F T F n H , 


til 9 

e 2 


e Q 

(9) 


FIG. 2: Profile of /ieg(r). The solid curve for r/ c (hl = 3), 
and the dashed one for r/b (hl =4). 


where n F = 3 (4) for eQ = e c (e&). As a bonus, Eq. Q) 
can be converted into the order-a 2 prediction to the 
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r] c —> 77 partial width, with the aid of the formula 
r [r] c —> 77 ] = E ^-\F(0)\ 2 M^ c . Our result is consistent 
with, but much more accurate than, the previous 0(a 2 ) 
result (2lj| . For completeness, we also include the “light- 
by-light” contribution previously neglected in [ill . 

It is widely known that perturbative expansion for 
F(0) exhibits poor convergence behavior. Taking /iA = 
to, the NNLO corrections overshoot the NLO correc¬ 
tions for both r) c and rjb, but the situation for the rj c 
is provokingly unsatisfying. The large negative prefac¬ 
tor of ln^A in 0 implies that lowering the NRQCD 
factorization scale would effectively dilute the NNLO 
correction. In facL in analyzing 0(a 2 ) correction to 
J/ip(T) — > e + e _ [20|, it was suggested that a better 
choice is ha = 1 GeV, the typical bound state scale for 
both charmonium and bottomonium. With this scale set¬ 
ting, r]b —> 77 receives the accidentally vanishing NNLO 
correction, but the convergence of perturbative expan¬ 
sion for r] c has not yet been improved to a satisfactory 
degree. 



Q 2 (GeV 2 ) 



FIG. 3: The transition form factor F(Q 2 ) normalized to F( 0 ) 
in variation with Q 2 : 7*770 (upper); 7 * 7 rji, (lower). The solid 
dots with error bars represent the BaBar measurements 0. 
The dotted curve stands for the LO NRQCD predictions, the 
(blue) dashed curve for the NLO predictions, while the solid 
curve for the NNLO predictions. The uncertainty bands as¬ 
sociated with the NNLO predictions are obtained by varying 
Hr from i( Q 2 + to 2 ) to 2(Q 2 + m 2 ). There are two distinct 
NNLO bands, the upper (green) one corresponds to ha = to, 
the lower (red) one to ha = 1 GeV. We suppress the (red) 
/x.R-uncertainty band for 75 since it is too narrow to be visible. 


In Fig. [3] we plot the ratio |F(Q 2 )/I ? (0)| against Q 2 
over a wide range, juxtaposing our predictions at various 
levels of accuracy in a s with the BaBar data [§]. For 
F(Q 2 ), we choose hr = \/Q 2 + m 2 in 0 by default; for 
-F(O) in 0, we fix hr = to. Using the latest PDG de¬ 
termination of the MS masses to c (to c ) = 1.28 GeV and 
TOfc(TOb) = 4.18 GeV j!5, we convert them into the 2-loop 
quark pole masses as m c = 1.68 GeV and mb = 4.78 GeV 
with the aid of the package Run Dec [4l|. We also use 
Run Dec to evaluate the running o^ s . To assess the un¬ 
certainty induced by the renormalization scale, we slide 
Hr from I/a/ 2 to \/2 times the default choice. For com¬ 
pleteness, we also include the “light-by-light” scattering 
contribution in our phenomenological analysis. 

As expected, the LO NRQCD prediction describes the 
measured 7 * 77 c form factor quite well. The agreement 
survives after adding the modest negative NLO pertur¬ 
bative correction. Since the NLO prediction is stable 
against the variation of hr, we suppress its uncertainty 
band in Fig. [3] 

It is counterintuitive that the NNLO prediction ap¬ 
pears to have stronger /independence with respect to the 
NLO result. This symptom can be attributed to the un¬ 
usually large |/^( r )l- The NNLO prediction turns out 
to be also sensitive to the NRQCD factorization scale 
Ha (This is mainly because we keep F(Q 2 )/F{ 0) un¬ 
expanded). To discern the uncertainty induced by ha, 
we include two sets of NNLO predictions in Fig. [3l one 
with ha = Tn (the ultimate UV cutoff scale of NRQCD) 
and ha = 1 GeV (roughly the scale of mv for charmo¬ 
nium and bottomonium). From the discussions following 
0 , we see that the former choice results in rather poor 
convergence for perturbative series. Indeed, it yields a 
prediction to |F(< 5 2 )/F( 0 )| that severely deviates from 
the BaBar data over the experimentally accessible Q 2 
range. For the latter scale setting, the discrepancy is 
considerably relieved, yet still pronounced, 5 a away from 
most data points. To reconcile the NRQCD prediction 
with the data, it is tempting to further decrease the value 
of ha, however one then faces the dilemma that NRQCD 
factorization would cease to make sense. 

In passing, it is worth mentioning that there exist 
different ways to “interpret” the NNLO corrections to 
F(Q 2 )/F(0), be., to expand this ratio strictly as power 
series in a s through the 2 nd order, or take the renor¬ 
malization scales associated with F(Q 2 ) and F( 0) in¬ 
terchangeably. In some scenarios, the NNLO correction 
turns out to be tiny and there reaches a perfect agreement 
between the “predictions” and data. A comprehensive 
survey reveals that, keeping F{Q 2 )/F{ 0) expanded or 
not is not the key that influences the agreement, rather, 
the key is rooted in whether one takes two hr s entering 
F(Q 2 ) and F( 0) equal or not. As a consequence of the 
deeply negative and flat (r), we find those predictions 
with equal- hr recipe will always yield decent agreement 
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with the data (no matter expanding F(Q 2 )/F(0 ) or not), 
those with unequal hr (like in our case) will generally 
result in severe discrepancy between the theory and the 
data. Note that, unlike the ratio of the S- wave quarko- 
nium decay rates as studied in [U, 24 1 , it is quite un¬ 
natural to take the hr in rj c production process equal to 
that in the r] c decay. They rather should be taken as the 
characteristic hard scale in respective production/decay 
processes. Moreover, since the experimental data are 
available separately for F(Q 2 ) and F(0), we can always 
convert the measured rj c —> 77 width to |U( 0 )|, then 
extract the NRQCD matrix element to the NNLO ac¬ 
curacy, then plug it to ([3]) to predict the F{Q 2 ) over a 
wide range of Q 2 . In computing the ratio F(Q 2 )/F(0), 
we would reproduce exactly the same plot as in Fig. [3] 
In our opinion, Fig. [3] remains the most honest way to 
present our NNLO prediction. It is important to note 
that the rise-then-drop shape of our prediction to F(Q 2 ) 
is qualitatively incompatible with the monopole shape 
observed experimentally. 

For 7*7 —> r)b, the size of NNLO correction to 
|F(Q 2 )/F( 0 )| seems to be modest even at p a = TOf,, and 
accidentally small for fi a = 1 GeV. Thus, the NRQCD 
factorization appears to work reasonably well for bot- 
tomonium. One has to await the next generation of high- 
energy e + e - collider to test our predictions. 

We tend to conclude that, within the reasonable ranges 
of ha and HRi the state-of-the-art NRQCD prediction 
utterly fails to account for the BaBar data on the 
7*7 Vc form factor. Since our predictions do not 

involve any adjustable nonperturbative parameters, and 
the disagreement prevails for a wide range of data points, 
this discrepancy poses a much more thorny challenge 
to theory than what was encountered in the e + e - —> 
J/ip + ri c 

A pressing question is how to resolve this disturbing 
puzzle within the confine of NRQCD. The relativistic 
correction to this transition form factor appears to be 
modest. Certainly it would be extremely illuminating to 
know the actual size of O(a^) correction, provided such 
a calculation were feasible in the foreseeable future. Fur¬ 
thermore, as was noted in 421, an outstanding short¬ 
coming of fixed-order NRQCD calculation for quarko- 
nium production is that, since the NRQCD matching 
coefficient contains several distinct scales Q , to, ha , it 
is a priori unclear how to optimally assign these scales 
to various a s . Without getting rid of this inherent scale 
ambiguity, it is difficult to make a sharp higher-order 
prediction out of the NRQCD framework. 

If the future Belle measurement of the 77 * —> r] c tran¬ 
sition form factor confirms the BaBar data, and if this 
discrepancy persists, serious doubt will be cast onto the 
usefulness of the NRQCD factorization approach when 
applied to charmonium exclusive production processes. 
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